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PREFACE 


This  report  documents  utility  code  developments  and  modifications  for  MODFLOW  and 
MODFLOWP  which  provide  new  model  set-up,  performance  monitoring  and  output  capabilities. 
Source  code  for  these  modifications  is  available  under  the  listing  for  this  report  through  the  home 
page  of  the  Center  for  Geophysical  Investigation  of  the  Shallow  Subsurface  at 
http://kihei.idbsu.edu/cgiss_pub.html.  This  report,  including  the  modified  MODFLOW  and 
MODFLOWP  codes,  are  public  domain  and  as  such  may  be  used  and  copied  freely.  The 
developed  programs  have  been  tested.  However,  no  warranty  is  given  that  the  code  is  completely 
error-free.  No  liability  is  assumed  for  any  damage  or  loss  that  may  result  from  use  of  our 
programs.  If  you  do  encounter  problems  with  the  code,  find  errors,  or  have  suggestions  for 
improvement,  please  contact  one  of  the  authors  at  the  addresses  below. 


Kangle  Huang 
Nuclear  Waste  Department* 
Earth  Science  Division 
MailStop  90-1116 
Lawrence  Berkeley  Laboratory 
Berkeley,  CA  94720 


Tel:  510-495-2454 
Fax:  510-486-5686 
Email:  khuang@lbl.gov 


TomClemo  Tel:  208-426-1416 

Warren  Barrash  Tel:  208-426-1229 

Center  for  Geophysical  Investigation  of 
the  Shallow  Subsurface 
Boise  State  University 
1910  University  Drive 
Boise,  ID  83725 


Email:  tomc@cgiss.idbsu.edu 
Email:  wb@cgiss.idbsu.edu 
Fax:  208-426-3888 


*  Present  address.  Kangle  Huang  created  the  MODFLOW  and  MODFLOWP  enhancements 
while  with  the  Center  for  Geophysical  Investigation  of  the  Shallow  Subsurface. 
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ABSTRACT 


This  report  documents  (1)  a  number  of  stand-alone  utilities  to  assist  with  input  to 
MODFLOW  and  (2)  code  modifications  internal  to  MODFLOW  and  MODFLOWP.  These  were 
developed  to  improve  ease  of  use  of  MODFLOW  and  MODFLOWP  for  modeling  groundwater 
flow  in  three-dimensional  heterogeneous  systems.  Several  utilities  to  assist  with  input  to 
MODFLOW  are  specifically  designed  for  use  with  the  Groundwater  Modeling  System  (GMS) 
pre-  and  post-processor  package. 

New  features  for  MODFLOW  include  routines  for:  ( 1 )  assignment  of  material  blocks  and 
material  types  to  these  blocks  within  a  given  three-dimensional  grid  mesh;  (2)  generation  of  the 
Block  Centered  Flow  (BCF)  file  for  heterogeneous  aquifers  with  material  properties  assigned  to 
blocks,  and  with  internal  calculation  of  leakance  and  thickness-factored  parameters  (T,  S);  (3) 
distribution  of  pumping  rate  for  individual  model  layers  for  a  well  that  draws  from  or  injects  into 
multiple  layers;  (4)  calculation  of  weighted  drawdown  for  an  observation  well  screened  over 
multiple  layers;  and  (5)  customizing  the  drawdown  output  file  to  contain  data  for  a  specified 
number  of  observation  wells,  rather  than  every  node  in  the  domain. 

For  MODFLOWP,  the  program  has  been  modified  to:  (1)  simplify  input  by  making  the 
.PAR  file  free  format;  (2)  allow  the  user  to  easily  change  the  number  and  combination  of 
parameters  for  estimation  in  a  given  run  by  modifying  one  line  rather  than  regenerating  the  .PAR 
file;  (3)  automatically  calculate  ROFFs,  COFFs,  and  TOFFs  for  observation  wells;  (4)  notify  the 
user  during  program  execution  when  lower  or  upper  parameter  ranges  are  exceeded;  and  (5) 
provide  monitoring  of  run  time. 
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1.  INTRODUCTION 


MODFLOW  (McDonald  and  Harbaugh,  1988)  is  perhaps  the  most  widely-used  numerical 
code  for  groundwater  modeling  due  to  its  ability  to  handle  complex  groundwater  flow  situations, 
including  three-dimensional  irregular  solution  domains,  heterogeneous  aquifers,  and  a 
comprehensive  range  of  boundary  conditions  and  flow  processes.  It  has  been  thoroughly  tested 
against  analytical  solutions  and  has  been  used  for  a  wide  variety  of  field  applications.  An  inversion 
code  based  on  MODFLOW,  MODFLOWP  (Hill,  1992),  provides  automated  estimation  of 
hydrologic  parameters.  We  have  been  using  both  MODFLOW  and  MODFLOWP  in  our 
modeling  of  pumping  tests  in  a  shallow  alluvial  aquifer  where  both  pumping  and  observation  wells 
are  partially  penetrating  wells  over  different  depth  intervals  in  the  aquifer.  Modeling  of  these  tests 
requires  a  three-dimensional  distribution  of  aquifer  parameters  to  achieve  matches  with  observed 
behavior  (Barrash  et  al.,  1995;  1997). 

This  report  is  an  outgrowth  of  these  efforts  to  model  heterogeneous  three-dimensional 
systems  and  presents  a  number  of  utilities  and  code  modifications  to  expedite  use  of  MODFLOW 
(Section  2)  and  MODFLOWP  (Section  3)  by  improving  input  and  output  ease  and  flexibility.  For 
MODFLOW,  some  of  these  utilities  are  designed  to  operate  in  conjunction  with  the  Department 
of  Defense  Groundwater  Modeling  System  (GMS),  a  comprehensive  pre-  and  post-processing 
software  package  for  MODFLOW,  as  well  as  other  common  groundwater  modeling  codes. 

In  particular,  this  report  presents  utilities  to  address  the  following  input  and  output  issues 
associated  with  three-dimensional  systems  and  pumping  test  applications  of  MODFLOW:  (1) 
assignment  of  material  blocks  and  material  types  to  these  blocks  within  a  given  grid  mesh;  (2) 
generation  of  the  Block  Centered  Flow  (BCF)  file  for  heterogeneous  aquifers  with  material 
properties  assigned  to  blocks  and  internal  calculation  of  leakance  and  thickness-factored 
parameters  (T,  S);  (3)  distribution  of  pumping  rate  for  individual  model  layers  for  a  well  that 
draws  from  or  injects  into  multiple  layers;  (4)  calculation  of  weighted  drawdown  for  an 
observation  well  screened  over  multiple  layers;  and  (5)  creation  of  a  drawdown  output  file  to 
contain  data  for  only  a  specified  number  of  observation  wells  rather  than  for  every  node  in  the 
model. 

For  MODFLOWP,  the  program  has  been  modified  to:  (1)  simplify  input  by  making  the 
.PAR  file  free  format;  (2)  allow  the  user  to  easily  change  the  number  and  combination  of 
parameters  for  estimation  in  a  given  run  by  modifying  one  line  rather  than  regenerating  the  whole 
.PAR  file;  (3)  automatically  calculate  ROFFs,  COFFs,  and  TOFFs;  (4)  notify  the  user  when 
parameter  values  go  beyond  upper  and  lower  limits  for  parameters  during  estimation  runs;  and  (5) 
monitor  run  time. 


2.  UTILITIES  FOR  USE  WITH  MODFLOW 
2.1  Automatic  Generation  of  Grid  Cells  and  Assignment  of  Cell  Properties 
2.1.1.  Problem 

Many  natural  geologic  materials  are  heterogeneous.  Accurately  characterizing  the  spatial 
distribution  of  hydraulic  properties  may  be  the  most  critical  step  in  modeling  groundwater  flow 


and  transport  in  strongly  heterogeneous  aquifers.  A  heterogeneous  aquifer  commonly  is 
characterized  by  dividing  the  aquifer  into  many  material  blocks  or  zones,  each  of  which  is 
assumed  to  be  homogeneous  (i.e.,  to  have  constant  hydraulic  properties  associated  with  the 
geologic  material  within  each  block  or  zone). 

GMS  defines  a  block  (3D)  or  zone  (2D)  by  selecting  a  group  of  cells,  and  it  then  assigns  a 
material  type,  MATID,  to  the  block  or  zone  according  to  user-provided  information.  (NOTE:  for 
brevity  in  this  report  we  will  henceforth  refer  only  to  blocks)  This  procedure  is  applied  to  define 
all  blocks  in  the  solution  domain.  Hydraulic  properties  of  each  grid  cell  are  then  identified 
through  the  MATID.  In  practice  this  procedure  can  be  very  time-consuming  and  error-prone, 
especially  when  the  model  structure  is  highly  heterogeneous. 

Not  only  must  three-dimensional  grids  be  generated  before  the  material  blocks  can  be 
defined,  but  commonly  the  model  is  modified  after  initial  runs  such  that  the  grid  must  be  changed, 
and/or  the  geometry  or  distribution  of  blocks  are  changed.  As  a  result,  previously  defined 
material  blocks  are  no  longer  useful  because  the  associated  cells  have  been  changed.  Hence,  the 
material  blocks  have  to  be  redefined  for  each  new  discretization.  Without  a  utility  to  automate 
this  process,  significant  amounts  of  time  are  needed  to  generate  each  new  or  modified  model 
structure,  introducing  the  potential  for  errors  to  be  incorporated. 

2.1.2,  Solution 

We  have  developed  the  program  BLOCK.FOR  to  accurately  define  the  material  blocks 
and  assign  the  MATIDs  for  all  cells.  Whenever  the  spatial  discretization  is  changed,  new 
MATIDs  can  be  obtained  by  running  the  BLOCK.FOR  program.  The  BLOCK.FOR  program  is  a 
stand-alone  utility  for  use  with  GMS  as  a  pre-processor  to  generate  the  three-dimensional  grid  file 
(filename. 3DG). 

BLOCK.FOR  works  in  the  following  way.  Assume  the  solution  domain  is  a  volume 
consisting  of  NLAY  model  layers,  each  of  which  has  NCOL  columns  and  NROW  rows. 
Coordinates  of  the  columns  and  rows  are  defined  by  X(z'),  i  =  l,...,NCOL+l  and  Y  (j),j  = 
l,...,NROW+l,  respectively  (Figure  1).  Also  assume  that  each  model  layer  (not  necessarily  the 
same  as  a  geologic  or  hydrostratigraphic  layer)  includes  several  types  of  materials,  and  each 
material  has  a  different  set  of  hydraulic  properties  that  is  constant  for  that  type  of  material.  On  a 
model-layer  by  model- layer  basis,  BLOCK.FOR  reads  information -on  the  number  and  locations  of 
blocks  and  the  material  types  associated  with  the  blocks.  To  reduce  the  input  effort,  the  user 
assigns  a  background  MATID  for  each  layer,  which  usually  would  be  the  material  type  that  is 
most  abundant  in  a  given  layer.  Thus,  the  user  need  only  define  the  specific  block  geometry  and 
MATIDs  for  the  remaining  portions  of  a  given  layer. 
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Figure  1.  Plan  view  of  material  blocks  described  in  the  example  application  (Section 
2.1.4).  Layer  1  is  at  the  top  and  Layer  2  is  at  the  bottom. 


Mat  1 

Mat  3 

Mat  2 

In  plan  view  each  block,  k,  in  a  given  layer  (i.e.,  k  =  1  to  NBLOCK,  or  total  number  of 
blocks  in  a  given  layer)  is  a  rectangular  area  that  is  uniquely  defined  by  the  coordinates  of  its 
lower  left  corner  (XLL(£),YLL(&))  and  upper  right  corner  (XUR(&),YTJR(/c)).  The  hydraulic 
properties  of  each  block  are  identified  by  its  material  type,  MATID(k).  With  this  information,  the 
aquifer  system  heterogeneity  can  be  built  into  the  model  using  BLOCK.FOR. 
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2.1.3.  Input  Instructions  for  the  BLOCK.FOR  program 

The  BLOCK.FOR  program  requires  two  input  files:  BLOCK.IN  and  a  .3DG  file.  These 
files  are  described  below  with  examples. 

2.1.3.1.  BLOCK.IN 

The  BLOCK.IN  file  defines  the  geometry  and  MATID  of  each  block  and  gives  the 
background  material  type  for  each  layer.  Thus,  there  are  NLAY  data  sets  in  BLOCK.IN,  one  for 
each  layer,  m  (i.e.,  m  =  1  to  NLAY).  Inputs  for  each  layer  are: 

1 .  MATID(m),  NBLOCK(m) 

Background  material  type  and  number  of  blocks  for  layer  m.  A  line  then  follows  for  each 
block,  k  (k  =  1  to  NBLOCK(m)),  in  layer  m: 

2.  MATID(fc),  XLL(Jt),  YLL(fc),  XUR(Jt),  YUR(fc) 

Materia]  type  and  coordinates  of  the  lower  left  (LL)  and  upper  right  (UR)  corners  of  block 

k.  If  this  layer  is  homogeneous  then  NBLOCK(m)  =  0.  Inactive  blocks  have  MATID(yt)  =  0.  It  is 
important  that  the  defined  blocks  in  a  layer  should  not  overlap,  although  they  will  have  corners 
and  sides  in  common  with  adjacent  blocks. 

2.1. 3.2.  .3DG  File 

The  .3DG  file  defines  the  solution  domain  and  its  discretization.  It  contains  the  following 
information: 

l .  NLAY,  NCOL,  NROW  Numbers  of  layers,  columns,  and  rows 

2.  X(i),  /=1,  NCOL+1  X  coordinates  of  nodes  in  a  layer 

3.  Y (j),j  =  1,  NROW+1  Y  coordinates  of  nodes  in  a  layer 

A  .3DG  file  can  be  obtained  directly  from  GMS  in  the  following  way:  when  a  solution 
domain  has  been  defined  in  GMS  (i.e.,  the  grid  has  been  formed),  save  the  grid  information  in  a 
file  using  the  GMS  command  sequence:  File— >  Save  and  provide  a  filename  to  go  with  the 
filespec  of  .3DG.  The  BLOCK.FOR  program  modifies  this  file  and  then  uses  the- same  filename 
for  the  output  file,  overwriting  the  original  file.  This  modified  new  filename. 3DG  file  has 
MATID  information  for  each  cell  -  which  is  what  is  needed  by  GMS  to  generate  the  model 
structure  (i.e.,  distribution  of  materials  in  a  given  three-dimensional  grid).  Then,  using 
BCF_IN.FOR  the  filename.3DG  file  can  be  used  to  generate  the  BCF  file  needed  by  MODFLOW 
(section  2.2).  Once  the  new  filename. 3DG  file  is  generated,  it  remains  valid  even  if  the  horizontal 
discretization  of  the  grid  is  changed  as  long  as  the  following  have  not  changed:  (1)  the  overall 
model  dimensions,  (2)  the  number  and  dimensions  of  layers,  and  (3)  the  locations  of  blocks. 

2.1.4.  Example  Application 

This  example  shows  how  the  BLOCK.FOR  program  works.  The  input  and  output  files 
for  the  example  are  also  illustrated.  First,  assume  a  region  of  an  aquifer  (surface  area  of  10  x  6 
length  units2)  is  discretized  into  blocks  bounded  in  the  x  direction  at  0,  2,  4,  7,  9,  and  10  length 
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units,  and  in  the  y  direction  at  0,  1,  3,  and  6  length  units.  The  model  of  the  aquifer  has  two 
layers:  the  upper  layer  is  5  units  thick  and  the  lower  layer  is  10  units  thick  (Figure  1).  Three 
materials  (MATID  1 , 2,  and  3)  are  identified  and  their  distributions  in  each  layer  are  shown  in 
Figure  1.  In  layer  1,  the  most  prevalent  material  (i.e.,  background  material)  is  MATID  1.  There 
are  two  blocks  defined  in  layer  1 :  one  consists  of  MATID2  and  is  defined  by  its  lower  left  corner 
(4,  1)  and  upper  right  corner  (9,  3);  the  other  block  consists  of  MATID3  and  is  defined  by  lower 
left  corner  (2,  0)  and  upper  right  corner  (7,  1).  In  layer  2,  the  background  material  is  MATID3 
and  there  are  two  material  blocks  as  in  layer  1 :  a  block  with  MATID  1  defined  by  lower  left 
corner  (0,  3)  and  upper  right  corner  (4,  6);  and  a  block  with  MATID2  defined  by  lower  left  corner 
(7,  0)  and  upper  right  corner  (10,  3). 

It  should  be  noted  that,  the  convention  used  here  for  locating  blocks  has  the  origin  in  the 
lower  left-hand  comer  of  a  layer  when  viewing  the  layer  in  plan  view.  However,  the  default  origin 
used  by  GMS  is  located  at  (100,  250,  -10)  for  the  example  of  Figure  1  and  Example  Files  1  and  2. 
It  should  be  noted  that  the  internal  GMS  coordinate  system  can  be  specified  by  the  modeler  such 
that  it  agrees  with  the  BLOCK.FOR  routine. 


Example  File  1.  Example  BLOCK.IN  file 


BLOCK.  IN 


Explanation 


1  2 

2  4.  1.  9.  3. 

3  2.  0.  7.  1. 
3  2 

1  0.  3.  4.  6. 

2  7.  0.  10.  3 


1=MATID1  for  background;  2=two  blocks  will 
2=MATID2,  lower  left  (4,1)  and  upper  right 
3=MATID3,  lower  left  (2,0)  and  upper  right 
3=MATID3  for  background,  2=two  blocks  will 
1=MATID1 ,  lower  left  (0,3)  and  upper  right 
2=MATID2,  lower  left  (7,0)  and  upper  right 


be  defined  in  following 

(9.3)  corners 
(7,1)  corners 

be  defined  in  following 
(4,6)  corners 

(10.3)  corners 


lines 


lines 


Below  is  an  example  .3DG  file  generated  when  a  model  grid  is  defined  in  GMS,  and  in 
which  MATIDs  are  not  yet  defined.  This  file  is  used  for  input  to  BLOCK.FOR. 


Example  File  2.  Example  .3DG  file 


Test.3DG 


GRID  three-dimensional 

TYPE  1 

IJK  -y  +x  -2 

ORIGIN  .  100000000000000E+03  . 250000000000000E+03  - . 100000000000000E+02 

ROTZ  O 

DIM  364 

•OOOOOOOOOOOOOOOE+OO 
•200000000000000E+01 
.400000000000000E+01 
.700000000000000E+01 
.900000000000000E+01 
. 100000000000000E+02 
.OOOOOOOOOOOOOOOE+OO 
•lOOOOOOOOOOOOOOE+Ol 
. 3 0000000000000 0E+01 
. 6000000000 0000 0E+01 
.000000000000000E+00 
■500000000000000E+01 
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. 100000000000000E+02 


The  output  file  is  Test. 3DG  which  has  been  overwritten  with  the  following  when 
BLOCK.FOR  is  run: 


Example  File  3.  Example  output  (overwritten)  .3DG  file 


TesUDG 


GRID  THREE-dimensional 
TYPE  1 

IJK  -y  +x  -z 

ORIGIN  .100000000000000E+03  . 250000000000000E+03  - . 100000000000000E+02 

ROTZ  0 

DIM  364 

.OOOOOOOOOOOOOOOE+OO 
.200000000000000E+01 
. 400000000000000E+01 
. 7 000000000 00000E+01 
. 900000000000000E+01 
. 100000000000000E+02 
.OOOOOOOOOOOOOOOE+OO 
.lOOOOOOOOOOOOOOE+Ol 
. 3 000000 00000000E+01 
.600000000000000E+01 
.000000000000000E+00 
. 500000000000000E+01 
.100000000000000E+02 

MAT 

111111122113311113333332233322 

This  output  .3DG  file  can  then  be  directly  read  into  GMS  to  completely  define  the  grid 
system  and  the  distribution  of  hydraulic  properties.  It  can  be  used  also  to  generate  a  portion  of 
the  MODFLOW  .BCF  file. 


2.2  Automatic  Generation  of  Block  Centered  Flow  (BCF)  Package 
2.2.1.  Problem 

In  modeling  three-dimensional  groundwater  flow  and  transport  in  heterogeneous  aquifers, 
the  values  of  hydraulic  properties  must  be  accurately  assigned  to  proper  grid  cells  in  the  numerical 
model.  This  task  can  be  time-consuming  and  prone  to  error  introduction,  even  with  the  aid  of  a 
powerful  pre-processor  such  as  GMS.  We  have  written  the  program  BCF_IN.FOR  to  efficiently 
prepare  the  Block  Centered  Flow  (BCF)  package  which  holds  all  the  information  describing  the 
flow  system  and  associated  hydraulic  properties.  The  values  of  some  hydraulic  properties  have  to 
be  calculated  (e.g.,  T,  S,  and  leakance)  and  then  input  to  a  selected  material  block  (group  of 
contiguous  cells  with  the  same  hydraulic  properties).  This  procedure  is  repeated  for  all  blocks. 
Completing  the  .BCF  input  file  manually  for  a  highly  heterogeneous  aquifer  system  takes  hours, 
or  even  a  day  or  more.  And  again,  the  laborious  process  must  be  repeated  to  generate  a  new 
.BCF  file  whenever  values  of  the  hydraulic  properties  or  block  dimensions  are  changed  during 
subsequent  model  runs. 
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2.2.2.  Solution 

The  basic  approach  of  BCF_IN.FOR  is  to:  (1)  directly  relate  the  geometric  information 
specifying  the  material  blocks  to  information  on  material  properties  (e.g.,  MATED,  K,  Ss  and  Sy); 
and  (2)  automatically  calculate  and  assign  values  of  hydraulic  properties  (T,  S,  leakance)  to  cells 
of  the  model.  Using  the  program  BCF_IN.FOR,  a  new  .BCF  file  can  be  generated  in  a  few 
seconds  when  the  values  of  hydraulic  properties  are  changed.  The  output  .BCF  file  from 
BCF_IN.FOR  can  be  directly  incorporated  into  the  MODFLOW  model  using  the  .BCF  input  file. 
BCFJN.FOR  is  an  independent  program.  However,  the  most  important  input  information  for 
this  program  comes  from  the  three-dimensional  grid  file  (filename. 3DG)  generated  by  GMS. 

2.2.3.  Input  Instructions  for  the  BCF_IN.FOR  program 

Two  input  files  in  addition  to  the  .BCF  file  are  required  by  the  BCF_IN.FOR  program:  a 
file  a  .3DG  file  and  a  file  called  MATERIAL.  IN 

2.2.3.1.  Three-Dimensional  Grid  File  (fllename.3DG) 

In  GMS,  after  three-dimensional  grids  have  been  generated  and  all  material  blocks  defined, 
the  three-dimensional  grid  file  can  be  obtained  by  choosing:  File  — >  Save  and  then  entering  the 
file  name. 

2.2.3.2.  MATERIAL.IN 

In  the  material  property  file  MATERIAL.IN,  list  each  type  of  material  using  one  line  for 
each.  In  each  line,  first  the  material  type  (MATID)  is  given  followed  by  K,  Sy  and  Ss  for  that 
material  type. 

2.2.4.  Example  Application 

We  demonstrate  using  the  BCF_IN.FOR  program  to  generate  the  .BCF  input  file  for  a 
hypothetical  heterogeneous  aquifer.  This  aquifer  is  composed  of  five  different  geological 
materials:  cobble-sand  mixture,  gravel,  clean  sand,  silty  sand,  and  silt.  These  materials  occur  in  a 
patchy  (rather  than  layered)  distribution  as  is  common  in  alluvial  aquifers.  The  simulated  volume 
is  about  100  m  x  100  m  in  area  and  40  m  thick.  The  solution  domain  is  discretized  into  14 
columns  and  14  rows  with  variable  grid  spacing  that  gradually  increases  from  the  center  of  the 
domain  (where  a  pumping  well  is  located)  to  the  boundaries.  The  vertical  dimension  is  divided 
into  three  model  layers  having  thicknesses  of  8.5  m,  18.5  m,  and  13  m  respectively,  from  bottom 
to  top  (Figure  2).  Material  blocks  are  defined  in  GMS  and  the  three-dimensional  grid  file  has 
been  saved  (select:  File— >  Save  and  give  the  filename  for  .3DG  file). 
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14  -  Rows 


Layers 


Cobbles  and  Sand 
Gravel 
Clean  Sand 
Silty  Sand 
Sandy  Silt 


B. 


Cobbles  and  Sand 
Gravel 
Clean  Sand 
Silty  Sand 
Sandy  Silt 


Figure  2.  Three-dimensional  view  and  plan  views  of  a  heterogeneous  aquifer  model  A.  Oblique 
view  of  a  three-layer  three-dimensionally  heterogeneous  aquifer  model  with  five  aquifer  material 
types.  B.  Plan  view  of  layer  1  (top  layer). 
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14  -< -  Rows 


C. 


A 


1  < 


Cobbles  and  Sand 
Gravel 
Clean  Sand 
Silty  Sand 
Sandy  Silt 


Columns 


>  14 


D. 


A 


w 


5 


o 

CE 


t 

■sT 


Cobbles  and  Sand 
Gravel 
Clean  Sand 
Silty  Sand 
Sandy  Silt 


1  < 


Columns 


>  14 


Figure  2.  Continued.  C.  Plan  view  of  layer  2  (middle  layer)  of  the  model.  D.  Plan  view  of  layer  3 
(bottom  layer)  of  the  model. 
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Following  is  the  .3DG  file  saved  by  GMS  for  this  example  application: 


Example  File  4.  Example.3DG  file 

EXAMPLE.  3  DG 


GRID3D 
TYPE  1 

IJK  -y  +x  -z 

ORIGIN  1 . 000000000000000e+002  2.000000000000000e+002  -4.000000000000000e+001 

ROTZ  0 

DIM  15  15  4 

0 . 000000000000000e+000 

1 . 000000000000000e+00 1 

2 . 000000000000000e+001 

2 . 800000000000000e+001 

3.600000000000000e+001 

4 . 100000000000000e+001 

4.600000000000000e+001 

4 . 900000000000000e+001 

5 . 400000000000000e+001 

5.900000000000000e+001 

6 . 7  00000000000000e+001 

7 . 500000000000000g+001 

8. 500000000000000e+001 

9 . 500000000000000e+001 

1 . 000000000000000e+002 

0 . 000000000000000e+000 

1 . 000000000000000e+001 

2 . 000000000000000e+001 

2 . 800000000000000e+001 

3 . 600000000000000e+001 

4. 100000000000000e+001 

4 . 600000000000000e+001 

4 . 900000000000000e+00 1 

5 . 400000000000000e+00 1 

5 . 900000000000000e+001 

6 . 700000000000000e+001 

7 . 500000000000000e+001 

8. 500000000000 000e+001  -----  —  - 

9 . 500000000000000e+001 

1 . 000000000000000e+002 

0 . 000000000000000e+000 

8. 500000000000000e+000 

2.700000000000000e+001 

4 . 000000000000000e+001 

MAT 

333331111113333333311111133333333111111133333331111112223 

33331111212223333311115555244411111555555444111111555 

514444111111111144441122212111444422222221114444222222211 

14333331111111144433333331111333331111155553333311111 

555533333111115555333331111111553333311111111133333111133 

33343333111133333433331111333334333311113333342555555 

555111425555555551114225555552211142222222222111422222222 

22111551111155552225511111222444455222224444444552222 

244444445522222444444444444444444444444444444444444444444 

44444441111111115555511111111155555111111111555551111 

11111555551111111115555511111111155555 
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The  first  part  of  this  input  file  consists  of  geometric  information  on  the  location  of  grid 
lines  for  rows,  columns,  and  layers.  The  second  part  contains  MATIDs  (in  this  case  for  five 
materials:  1  through  5),  one  for  each  cell  in  the  model  grid  and  each  corresponding  to  the  MATID 
as  specified  by  the  example  .MAT  file  described  below. 


Example  File  5,  Example  .MAT  file 

EXAMPLE  MAT 


1  0.02  0.2  5.0e-04 

2  0.125  0.21  7.0e-04 

3  0.085  0.18  6.0e-05 

4  0.01  0.15  9.0e-05 

5  0.01  0.23  2.0e-04 


The  original  .BCF  file  that  will  be  taken  as  input  by  BCF_IN.FOR  also  is  obtained  from 
GMS  when  the  MODFLOW  model  domain  is  generated  in  GMS  and  saved  there.  The  data  in  the 
original  .BCF  file  will  just  be  read  in,  and  then  this  file  will  be  overwritten  as  the  new  .BCF  file 
with  the  same  name,  which  now  serves  as  the  output  file.  The  following  is  the  initial  .BCF  file 
(input  to  BCF_IN.FOR): 


Example  File  6.  Example  .BCF  input  file 

EXAMPLE  .BCF 


0  0  -888.0000 

0  1.0000000 

1 

0 

u 

0  1.0000000 

11  1.0000000 

(10G15.6) 

-1 

10.000000  10.000000 

8.000000 

8.000000 

5.000000 

5.000000 

3.000000 

5.000000 

5.000000  8.000000 

8.000000  10.000000 

10.000000 

5.000000 

11  1.0000000 

(10G15.6) 

-1 

5.000000  10.000000 

10.000000 

8.000000 

8.000000 

5. 000000 

5.000000 

3.000000 

5.000000  5.000000 

8.000000  8.000000 

10.000000 

10.000000 

0  0.0 

0  0.0 

0  0.0 

0  0.0 

0  0.0 

0  0.0 

0  0.0 

0  0.0 

0  0.0 


Shown  below  is  the  output  .BCF  file  (note  it  has  the  same  name  as  the  input  file) 
generated  by  running  the  BCF_IN.FOR  program. 


Example  File  7.  Example  .BCF  output  file 

EXAMPLE  BCF 


0  0  -888.0000  0  1.0000000 


1  0  0 

0 

1.0000000 

0 

11 

1.0000000 

(10G15.6) 

-1 

10.0000 

10.0000 

8.00000 

8.00000 

5.00000 

8.00000 

8.00000 

10.0000 

10.0000 

5.00000 

11 

1.0000000 

(10G15.6) 

-1 

5.00000 

10.0000 

10.0000 

8.00000 

5.00000 

5.00000 

8.00000 

8.00000 

10.0000 

10.0000 

11 

1.0000 

(10G11.4) 

0 

.1800 

.1800 

.1800 

.1800 

.1800 

.2000 

.1800 

.1800 

.1800 

.1800 

.2000 

.2000 

.2000 

.2000 

.2000 

.1800 

.1800 

.1800 

.2000 

.2000 

.1800 

.1800 

.1800 

.1800 

.1800 

.2000 

.2000 

.2000 

.2100 

.2100 

.1800 

.2000 

.2000 

.2000 

.2000 

.1800 

.1800 

.1800 

.1800 

.1800 

.2300 

.2300 

.2300 

.2100 

.1500 

.2000 

.2000 

.2300 

.2300 

.2300 

.1500 

.2000 

.2000 

.2000 

.2000 

.2300 

.2000 

.1500 

.1500 

.1500 

.2000 

.2000 

.2000 

.2000 

.2000 

.2000 

.2000 

.2100 

.2100 

.2100 

.1500 

.1500 

.1500 

.1500 

.2100 

.2100 

.2000 

.2000 

.2000 

.1500 

.2100 

.2100 

.2100 

.2100 

.2100 

.1800 

.1800 

.1800 

.1800 

.2000 

.2000 

.2000 

.1500 

.1500 

.1500 

.1800 

.1800 

.2000 

.2000 

.2000 

11 

1.0000 

(10G11.4) 

0 

.8500E-01 

.8500E-01 

.8500E-01 

.8500E-01 

.8500E-01 

•2000E-01 

.8500E-01 

.8500E-01 

.8500E-01 

.8500E-01 

.2000E-01 

.2000E-01 

.2000E-01 

.2000E-01 

.2000E-01 

.8500E-01 

.8500E-01 

•8500E-01 

.2000E-01 

.2000E-01 

.8500E-01 

.8500E-01 

.8500E-01 

.8500E-01 

.8500E-01 

•2000E-01 

.2000E-01 

.2000E-01 

.1250 

.1250 

.8500E-01 

.2000E-01 

.2000E-01 

.2000E-01 

.2000E-01 

.8500E-01 

.8500E-01 

.8500E-01 

.8500E-01 

.8500E-01 

. 1000E-01 

.1000E-01 

.lOOOE-Ol 

.1250 

.  1000E-01 

.2000E-01 

.2000E-01 

. 1000E-01 

. 1000E-01 

.1000E-01 

. 1000E-01 

.2000E-01 

.2000E-01 

.2000E-01 

.2000E-01 

. 1000E-01 

.2000E-01 

. 1000E-01 

. 1000E-01 

. 1000E-01 

.2000E-01 

.2000E-01 

.2000E-01 

.2000E-01 

.2000E-01 

.2000E-01 

.2000E-01 

.1250 

.1250 

.1250 

.1000E-01 

. 1000E-01 

.  1000E-01 

.1000E-01 

.1250 

.1250 

.2000E-01 

.2000E-01 

.2000E-01 

. 1000E-01 

.1250 

.1250 

.1250 

.1250 

.1250 

.8500E-01 

•8500E-01 

.8500E-01 

.8500E-01 

•2000E-01 

.2000E-01 

.2000E-01 

. 1000E-01 

.  1000E-01 

.  1000E-01 

.8500E-01 

.8500E-01 

.2000E-01 

.2000E-01 

.2000E-01 

0 

-13.0000 

11 

1.0000 

(10G11.4) 

0 

.5397E-02 

.5397E-02 

.5397E-02 

.5397E-02 

.5397E-02 

.8000E-03 

.9985E-03 

.9985E-03 

.9985E-03 

.5397E-02 

1  0 


5.00000 

5.00000 

3.00000 

5.00000 

.000000 

8.00000 

5.00000 

5.00000 

3.00000 

.000000 

Sy  of  unconfined  layer  1 

.2000 

.2000 

.2000 

.2000 

.2000 

.1800 

.1800 

.1800 

.1800 

.2000 

.1800 

.1800 

.1800 

.1800 

.1800 

.2000 

.2000 

.2000 

.2000 

.2000 

.1800 

.1800 

.2000 

.2000 

.2000 

.2100 

.1800 

.1800 

.1800 

.1800 

.2100 

.2000 

.2100 

.2100 

.2100 

.2000 

.2000 

.2000 

.2000 

.2300 

.1500 

.1500 

.2000 

.2000 

.2000 

.2300 

.2300 

.2300 

.1500 

.1500 

.2000 

.2000 

.2300 

.2300 

.2300 

.1500 

.2000 

.2000 

.2000 

.2000 

.2000 

.1500 

.1500 

.1500 

.1500 

.2000 

.2100 

.2000 

.2000 

.2000 

.2100 

.2100 

.2100 

.2100 

.2100 

.1500 

.1500 

.1500 

.2100 

.2100 

.2000 

.2000 

.2000 

.1500 

.1800 

.2000 

.2000 

.2000 

.2000 

.2000 

.1800 

.1800 

.1800 

.1800 

.1800 

.2000 

K  of  unconfined  layer  1 

.2000E-01 

.2000E-01 

.2000E-01 

.2000E-01 

.2000E-01 

.8500E-01 

.8500E-01 

.8500E-01 

.8500E-01 

.2000E-01 

.8500E-01 

.8500E-01 

.8500E-01 

.8500E-01 

•8500E-01 

.2000E-01 

.2000E-01 

.2000E-01 

.2000E-01 

.2000E-01 

.8500E-01 

.8500E-01 

.2000E-01_. 

.  .2000E-01 

.2000E-01 

.1250 

.8500E-0i 

.8500E-01 

T8500E-01 

.8500E-01 

.1250 

•2000E-01 

.1250 

.1250 

.1250 

.2000E-01 

.2000E-01 

.2000E-01 

.2000E-01 

.lOOOE-Ol 

.  1000E-01 

.1000E-01 

.2000E-01 

.2000E-01 

.2000E-01 

.1000E-01 

.1000E-01 

. 1000E-01 

.  1000E-01 

. 1000E-01 

.2000E-01 

.2000E-01 

•1000E-01 

.1000E-01 

. 1000E-01 

.  1000E-01 

.2000E-01 

.2000E-01 

.2000E-01 

■2000E-01 

.2000E-01 

. 1000E-01 

.  1000E-01 

.  1000E-01 

. 1000E-01 

.2000E-01 

.1250 

.2000E-01 

•2000E-01 

.2000E-01 

.1250 

.1250 

.1250 

.1250 

.1250 

.1000E-01 

.  1000E-01 

.  1000E-01 

.1250 

.1250 

.2000E-01 

.2000E-01 

.2000E-01 

.  1000E-01 

.8500E-01 

.2000E-01 

.2000E-01 

.2000E-01 

.2000E-01 

.2000E-01 

.8500E-01 

.2000E-01 

•8500E-01 

.8500E-01 

.8500E-01 

.8500E-01 

Bottom  elevation  of  unconfined  layer  1 
Leakance  of  unconfined  layer  1 


.1270E-02  . 1270E-02  .1270E-02  .1270E-02  .1270E-02 
.5397E-02  .5397E-02  .5397E-02  .5397E-02  .1270E-02 
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.  1270E-02 

.  1270E-02 

. 1270E-02 

•1270E-02 

■8000E-03 

.9985E-03 

.9985E-03 

.9985E-03 

.5397E-02 

. 5397E-02 

. 5  397  E - 02 

.5397E-02 

.5397E-02 

. 1270E-02 

.  1270E-02 

. 1270E-02 

.  1270E-02 

. 1270E-02 

•8000E-03 

.8000E-03 

.9985E-03 

.9985E-03 

.5397E-02 

.5397E-02 

.5397E-02 

.5397E-02 

.5397E-02 

. 1270E-02 

.  1270E-02 

. 1270E-02 

. 1270E-02 

.  1270E-02 

. 1270E-02 

. 1944E-02 

.1024E-02 

. 1024E-02 

.5397E-02 

.5397E-02 

•5397E-02 

.5397E-02 

.5397E-02 

.  1270E-02 

. 1270E-02 

. 1270E-02 

. 1270E-02 

.1944E-02 

.  1270E-02 

. 1944E-02 

.  1944E-02 

. 1944E-02 

■5397E-02 

.5397E-02 

.5397E-02 

.5397E-02 

.5397E-02 

. 1270E-02 

.  1270E-02 

. 1270E-02 

.  1270E-02 

. 1318E-02 

. 131 8E - 02 

•1318E-02 

. 1318E-02 

.6218E-02 

•6349E-03 

. 1318E-02 

.  1318E-02 

.2305E-02 

.2305E-02 

. 1270E-02 

. 1270E-02 

. 1270E-02 

.8989E-03 

. 1318E-02 

.  1318E-02 

..1318E-02 
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2.3.  Automatic  Generation  of  the  Well  Package  for  Wells  Open  to  Multiple  Layers 

2.3.1.  Problem 

A  pumping  well  may  be  screened  across  several  geologic  or  hydrostratigraphic  layers  having 
different  hydraulic  properties,  or  several  model  layers  in  a  numerical  model  such  as  MODFLOW  that 
have  the  same  or  different  hydraulic  properties.  The  pumping  rate  of  a  well  is  comprised  of  the  sum  of 
the  flow  rates  of  all  the  geologic  or  model  layers  contributing  to  the  well  discharge.  In  MODFLOW 
the  pumping  rate  (the  strength  of  the  sink  term)  must  be  specified  for  each  individual  model  layer.  It  is 
apparent  that  the  pumping  rate  for  a  well  in  an  individual  layer  is  a  fraction  of  the  total  pumping  rate  of 
the  well.  In  practice,  fractional  pumping  rates  for  individual  layers  have  been  approximated  (e.g., 
Javandel  and  Witherspoon,  1969;  Molz  et  al.,  1990;  McDonald  and  Harbaugh,  1998  eq.  68)  by 
apportioning  flow  based  on  the  relative  transmissivity  of  the  contributing  layers  (Figure  3).  Currently, 
the  MODFLOW  user  must  calculate  the  individual  pumping  rates  per  well  cell,  representing  the 
fractional  contribution  for  a  given  layer,  for  each  time  period  that  has  a  different  pumping  rate.  The 
calculations  must  be  repeated  if  model  parameter  values  or  grid  configuration  are  changed. 

Performing  such  calculations  by  hand  can  be  time-consuming  and  error-prone. 

2.3.2.  Solution 

A  simple  algorithm  is  implemented  in  the  program  WELL_Q.FOR  to  apportion  flow  rates  by 
layer  based  on:  (1)  values  of  T  designated  for  active  pumping  cells;  and  (2)  total  Q  for  a  given  period. 
The  assumption  behind  this  approach  is  that  the  drawdown  in  the  pumping  well  is  uniform  throughout 
the  well  so  that  an  aquifer  or  model  layer  having  a  relatively  large  transmissivity  value  contributes  a 
higher  flow  rate  to  the  total  pumping  rate  than  layers  with  lower  T.  Based  on  this  assumption,  the 
fraction  of  the  total  pumping  rate  contributed  by  a  well  in  an  individual  model  layer  can  be  calculated 
by  the  following  equation: 


n 

ETj 

J=  i 


where  n  is  the  number  of  model  layers  contributing  to  the  well  (for  example,  n  -  5  for  the  well 
shown  in  Figure  3),  i  is  the  layer  for  which  Qj  of  the  well  is  being  calculated  and  the  summation,  Tj; 
gives  the  combined  transmissivity  of  all  layers  open  to  the  well. 

2.3.3.  Input  Instructions  for  the  WELL_Q.FOR  Program 

Two  input  files  must  be  provided  for  the  WELL_Q.FOR  program:  WELL.IN  and  Q_TIME. 

2.3.3.I.  WELL.IN 

This  file  (see  Example  File  8)  provides  the  model  layer  number  having  an  active  well  cell,  the 
hydraulic  conductivity  (K)  of  the  layer,  and  the  thickness  of  this  layer.  For  each  multilayer  well,  the 
WELL.IN  file  has  one  line  for  each  layer  in  the  model  that  has  an  active  well  cell. 
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2.3.3.2.  Q_TIME 

In  this  file  (see  Example  File  9),  NCOL  (column  number)  and  NROW  (row  number)  of  the 
pumping  well  are  given  on  the  first  line.  Then,  one  line  is  given  for  each  stress  period  with  the 
first  number  on  the  line  being  the  stress  period  sequence  number,  and  the  second  number  on  the 
line  being  the  pumping  rate  for  that  stress  period  in  consistent  units.  The  pumping  rate  value  in 
Q_TIME  follows  the  MODFLOW  sign  convention  with  negative  value  for  pumping  (sink),  or 
positive  values  for  injection  (source). 

2.3.4.  Example  Application 

This  example  demonstrates  the  use  of  the  WELL_Q.FOR  program  for  generating  the  well 
package  for  a  MODFLOW  model  with  a  well  pumping  in  a  shallow  heterogeneous  aquifer.  A 
model  domain  of  250  ft  x  250  ft  x  58  ft  was  discretized  into  65  rows,  65  columns,  and  28  layers. 
The  pumping  well  consists  of  cells  in  three  model  layers  (layers  21  to  23)  in  NCOL  =  33  and 
NROW  =  33.  As  in  many  real  circumstances,  the  pumping  rate  in  this  model  scenario  fluctuates 
at  early  time.  To  account  for  this  early-time  variable  pumping  rate,  six  stress  periods  are  used. 
Input  files  for  this  example  are  given  below: 

Ja 


■  Cobbles  and  Sand 
Gravel 

H'  Clean  Sand 
Silty  Sand 
Sandy  Silt 

Figure  3.  Schematic  diagram  of  a  pumping  well  screened  on  five  layers 


Q1,T1 


Q2.T2 
Q3,  T3 

04,  T4 

05,  T5 


Example  File  8.  Example  WELL.IN  file 
WELL.  IN 

21  0.125  1.25 
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22  0.225  2.0 

23  0.225  1.0 

Example  File  9.  Example  Q_TIME  file 
O  TIME 

33  33 

1  -0.35 

2  -1.200000 

3  -1.800000 

4  -2.47000000 

5  -0.200000 

6  -2.4700000 


The  output  file  of  the  WELL_Q.FOR  program  is  the  well  package  file  for  MODFLOW: 


Example  File  10.  Example  VARIABLE. WEL  file 


VARIABLE.  WEL 


3  0 

3 


21 

33 

33 

-.065789 

22 

33 

33 

-.189474 

23 

33 

33 

-.094737 

3 

21 

33 

33 

-.225564 

22 

33 

33 

-.649624 

23 

33 

33 

-.324812 

3 

21 

33 

33 

-.338346 

22 

33 

33 

-.974436 

23 

■a 

33 

33 

-.487218 

21 

33 

33 

-.464286 

22 

33 

33 

-1.337143 

23 

33 

33 

-.668571 

3 

21 

33 

33 

-.037594 

22 

33 

33 

-.108271 

23 

33 

33 

-.054135 

3 

21 

33 

33 

-.464286 

22 

33 

33 

-1.337143 

23 

33 

33 

-.668571 

2.4.  Output  of  Drawdown  for  Observation  Wells 
2.4.1.  Problem 

MODFLOW  output  provides  drawdown  values  for  every  active  node  in  the  model  domain  at 
each  time  step  of  a  model  run.  Although  GMS  has  been  designed  to  be  able  to  retrieve  drawdown 
values  for  an  observation  point  from  MODFLOW  model  runs  through  the  GAGE  function,  it  is 
required  that  model-calculated  drawdown  be  saved  for  every  time  step.  As  a  result,  the  saved 
drawdown  output  file,  even  in  binary  form,  can  be  quite  large;  ASCII  Text  files  on  the  order  of  50  to 
100  MB  are  not  unusual  for  a  three-dimensional  model  having  a  hundred  thousand  cells  or  more.  For 


17 


some  modeling  activities  such  as  matching  drawdown  behavior  at  pumping  and/or  observation  wells, 
simulated  drawdown  results  are  needed  at  these  particular  wells  rather  than  throughout  the  full  flow 
domain  -  so  the  complete  set  of  drawdown  results  for  a  given  run  would  not  be  needed. 

2.4.2;  Solution 

To  reduce  the  amount  of  writing  to  disk  and  the  size  of  the  drawdown  file,  we  developed  the 
OBSERV.FOR  routine  for  recording  the  drawdown  versus  time  data  only  at  specified  observation 
points  or  cells,  rather  than  at  all  cells.  This  is  implemented  within  MODFLOW,  and  is  activated 
through  minor  modifications  to  the  GMS  super  file  and  the  MODFLOW  .BAS  file  (described  in  detail 
below  in  section  2.4.6.).  The  user  specifies  the  needed  observation  points  using  either  row,  column, 
and  layer  specifications  or  absolute  coordinates.  The  OBSERV.FOR  routine  uses  the  coordinate 
transformation  specification  contained  in  the  MODFLOW  .MFS  file  to  convert  from  coordinate 
positions  to  row,  column,  layer  indices. 

2.4.3  Input  Instructions  for  OBSERV.FOR 

2.4.3.1  The  .OBI  File 

The  inputs  for  the  OBSERV.FOR  routine  are  the  locations  of  well  cells  of  interest  placed  in  a 
file  commonly  using  the  suffix  .OBI.  The  locations  of  well  cells  can  be  represented  by  either  the  cell 
row,  column,  layer  indices  or  the  x,  y,  z  coordinates  of  an  observation  point,  or  well  cells.  Wells  that 
intercept  multiple  model  layers  have  cell  location  data  for  each  layer.  The  format  for  the  .OBI  file  is: 

NOBS,  OPTION 
For  each  NOBS  well  layer  if  OPTION  =  0 
I  J  K 

For  each  NOBS  well  if  OPTION  <  0 
X  Y  Z 

For  each  NOBS  well  if  OPTION  >  0 
X  Y  K 


Where  _ 

NOBS  -  Number  of  observation  locations 
OPTION  -  0  ->  Integer  row,  column  and  layer  specification 

<0  ->  x  ,  y,  z  coordinate  positions 

>0  ->  x,  y  coordinate  locations  with  K  as  an  integer  layer  specification 
I,  J,  K  INTEGER  row,  column  and  layer  specification  respectively 

X,  Y,  Z  REAL  coordinates 

2.4.4.  Application  Example 

The  use  of  the  OBSERV.FOR  routine  is  illustrated  using  test  case  1  of  the  MODFLOWP 
user’s  manual  (Hill,  1992).  The  pumping  well  is  located  at  row  9  column  10.  The  domain  contains 
two  layers  with  18  rows  and  columns.  The  horizontal  cell  widths  are  1,000  ft  on  each  side.  The  layers 
are  50  ft  thick  separated  by  a  10  ft  thick  aquitard  that  is  not  explicitly  represented  in  the  model.  The 
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well  penetrates  both  model  layers.  The  default  GMS  axis  definitions  of  x  corresponding  to  J  (columns) 
and  y  corresponding  to  I  (rows)  are  used.  X  increases  with  increasing  J  and  y  decreases  with 
increasing  I.  These  relationships  of  the  coordinate  axes  to  the  row  and  column  indices  are  defined  in 
the  GMS  .MFS  MODFLOW  superfile  (See  Example  file  15).  Below  is  the  .OBI  input  file  for 
OBSERV.FOR  using  cell  numbers  for  this  example: 

Example  File  11.  Example  .OBI  file  using  cell  numbers 

2.  0 
9.  10.  1 
9.  10,  1 


And  below  is  the  corresponding  .OBI  file  for  the  same  well  cells  if  the  coordinates  of  the  well 
cells  are  used  instead  of  the  cell  numbers: 

Example  File  12.  Example  .OBI  file  using  coordinates 

2.-1 

9500.,  9500.  ,  -40. 

9500.,  9500.  ,  -75. 


All  inputs  in  the  .OBI  file  are  free  format.  If  coordinates  are  used,  they  need  only  lie  within  the 
cell;  they  do  not  have  to  exactly  match  the  coordinates  of  the  cell  node.  The  well  coordinates  option 
may  be  most  useful  when  the  model  grid  geometry  is  irregular  and/or  when  the  geometry  is  changed 
for  subsequent  model  runs  -  because  then  one  does  not  need  to  update  the  row,  column,  and  layer 
numbers  for  the  observation  points. 

The  output  using  Example  file  12  is  listed  below  as  Example  file  13.  The  input  information  and 
GMS  coordinate  transformation  are  listed  at  the  top.  The  first  column  well  cell  data  is  the  time  step 
(KSTP),  second  the  stress  period  (KPER),  third  the  elapsed  time  since  the  start  of  pumping 
(T_TOTAL).  The  next  columns  are  calculated  drawdown  for  each  well  cell.  The  cell  numbers  are 
listed  on  the  top  of  each  column. 


2.4.5.  Weighted  Drawdown  Based  on  Transmissivity 

Because  the  model  equivalent  of  field-observed  drawdown  is  a  combined  drawdown  of  the 
cells  of  all  the  layers  open  to  the  well,  a  weighted  drawdown  has  to  be  calculated  for  comparison  of 
model  and  observed  results.  This  is  done  by  calculating  a  weighted  drawdown  for  the  pumping  or 
observation  well  based  on  the  values  of  transmissivity  (Tp  of  the  involved  cells  (e.g.,  Javandel  and 
Witherspoon,  1969).  For  example,  if  a  well  is  screened  in  n  cells,  the  drawdown  in  that  well  (s)  may 
be  calculated  as  the  sum  of  the  weighted  drawdown  in  the  contributing  n  cells: 


s 


E 

j=  i 


SJTJ 


where  T  is  the  sum  of  the  transmissivity  (Tj)  values  of  all  screened  cells.  Note  that  the  cells  are  not 
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necessarily  vertically  adjacent  to  each  other  -  the  layers  open  to  a  well  can  be  separated  by  blank  or 
cased  intervals.  However,  they  should  have  the  same  horizontal  coordinates  (i.e.,  their  column  and 
row  numbers,  respectively,  should  be  the  same). 


Example  File  13.  Example  of  .OBO  file 


Observation  location  drawdown 
2  cells  using  option  -1 
9500.  9500.  -40.00 

9500.  9500.  -70.00 

IJK  -y  +x  -z 

0RI6  O.OOOOOOE+OO  0.000000E+00  -100.000000 

R0TZ  0.000000E+00 

Observation  cell  drawdown 


KSTP 

KPER 

T_T0TAL 

R0W=  9  R0W=  9 

C0L=  10  C0L=  10 

LAY=  1  LAY=  2 

1 

1 

1.000 

■9853E-03 

•9992E-02 

1 

2 

4.000 

. 3942E-02 

•3991E-01 

1 

3 

10.00 

•9859E-02 

•9953E-01 

1 

4 

23.12 

■2281E-01 

.2288 

2 

4 

38.86 

. 3838E-01 

.3824 

3 

4 

57.74 

.5709E-01 

.5645 

4 

4 

80.41 

■7959E-01 

.7800 

5 

4 

107.6 

.1067 

1.034 

6 

4 

140.2 

.1392 

1.333 

7 

4 

179.4 

.1784 

1.684 

8 

4 

226.4 

.2256 

2.093 

9 

4 

282.8 

.2825 

2.567 

Weighted  average  drawdown  of  multilayer  wells 

KSTP 

KPER 

T_T0TAL 

R0W=  9 
C0L=  10 

1 

1 

1.000 

•4181E-02 

1 

2 

4.000 

. 1671E-01 

1 

3 

10.00 

■4168E-01 

1 

4 

23.12 

■9590E-01 

2 

4 

38.86 

.1604 

3 

4 

57.74 

.2371 

4 

4 

80.41 

.3281 

5 

4 

107.6 

.4358 

6 

4 

140.2 

.5630 

7 

4 

179.4 

.7126 

8 

4 

226.4 

.8881 

9 

4 

282.8 

1.093 

2.4.6.  Notes  on  Usage  with  GMS 

The  OBSERV.FOR  routine  is  included  in  the  modified  version  of  MODFLOW  that 
accompanies  this  report.  To  activate  this  modified  MODFLOW  code,  the  following  additional  files 
and  modifications  must  be  made: 


1.  Establish  the  input  file  named:  prefix. OBI  (see  section  2.4.4  above) 

2.  Edit  the  MODFLOW  super  file,  prefix. MFS,  generated  by  GMS  in  the  following  ways: 
a)  Add  the  following  two  lines  in  the  file  list: 

OBIN  31  prefix. OBI 
OBOT  32  prefix.OBO 

NOTE:  In  these  two  new  lines,  31  and  32  are  FORTRAN  unit  numbers;  any 
numbers  can  be  used  as  long  as  they  do  not  conflict  with  other  assigned  unit 
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numbers.  Besides  those  listed  in  the  .BAS  file  and  .MFS  files,  FORTRAN  unit 
numbers  98  and  99  are  reserved  for  a  temporary  scratch  file  and  the  .MFS  file 
respectively. 

b)  Put  single  quotation  marks  around  the  word:  ORIG  (see  Example  File  14  below). 

3.  Edit  line  4  of  the  .BAS  file  by  activating  IUNIT  number  17  for  the  input  file  (31)  and  IUNIT 
number  21  for  the  output  file  (32)  to  be  used  by  .OBI  and  .OBO  files,  respectively  (see  Example  File 
15  below). 

The  following  example  is  the  part  of  the  GMS  super  file  for  MODFLOW  which  illustrates  the 
changes  described  above: 

Example  File  14.  Example  modified  .MFS  file 

MODSUP 

IJK  -y  +x  -z 

LIST  26  "two_layer_trans.out" 

BAS1  1  "two_layer_trans.bas" 

BCF3  11  "two_layer_trans.bcf" 

0UT1  10  "two_layer_trans.oc" 

HEAD  -30  "two_layer_trans.hed" 

DRAW  -35  "two_l ayertrans .drw" 

CCF  -40  "two_layer_trans40.ccf" 

PCS2  12  "twolayertrans .peg" 

RIV1  15  "two_layer_trans.riv“ 

WEL1  13  "two_layer_trans .wel " 

GHB1  16  "two_layer_trans.ghb" 

RCH1  20  “two_layer_trans.rch" 

MT  three-dimensional  -29  "two_layer_trans.hff" 

OBIN  31  “ two_l ay er_tran  s .obi" 

0B0T  32  “two_layer_trans.obo" 

'ORIG'  0 . 000000000000000e+000  0.000000000000000e+000  -1.00000000000000e+002 
ROTZ  0 

LAYER  0  50.000000  (10G15.6) 

DMAT  1 
MPSC  10  0 
MTRN  1  0  0 
MHC  1  0  0 
MBE  1  0  0 
MLK  1  0  0 

MSSC  10  0  _ _  _ . 

MTE  1  0  0 
MWDF  10  0 


Example  File  15.  Top  of  modified  .BAS  file 

MODULAR  MODEL  -  TW0_LAYER  EXAMPLE  PROBLEM 

2  18  18  4  4 

u  13  o  15  o  o  16  20  o  o  o  io  12  o  o  o  3i  o  o  o  32  29  o  o  @  @  @  Modified  line 


0 

0 

0 


1 

1 

1 


-999.0000 


1  1.0000000 


(10G15.6) 


0 
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2.4.7.  Notes  on  Usage  with  MODFLOWP 

At  the  time  this  report  was  written,  GMS  used  the  9/1/87  version  of  MODFLOW  88  and  the 
latest  version  of  MODFLOWP  (3.2)  used  the  5/23/96  version  of  MODFLOW  96.  One  of  the 
differences  between  MODFLOW  88  and  MODFLOW  96  is  the  method  used  to  open  files.  To  activate 
the  OBSERV.FOR  routine  in  the  modified  MODFLOWP  program  do  the  following: 

1 .  Establish  the  input  file  named:  prefix.OBl  (see  section  2.4.4  above) 

2.  Add  the  following  two  lines  in  the  NAME  FILE: 

OB  IN  31  p at  h\p refix.  OBI 
OBOT  32  pathy?re/«c.OBO 

NOTE:  In  these  two  new  lines,  31  and  32  are  FORTRAN  unit  numbers;  any 
numbers  can  be  used  as  long  as  they  do  not  conflict  with  other  assigned  unit 
numbers.  Besides  the  FORTRAN  unit  numbers  listed  in  the  NAME  FILE,  unit 
numbers  0  and  96  through  99  are  reserved  in  MODFLOWP. 


3.  MODIFICATIONS  TO  MODFLOWP 

The  above  sections  present  our  modifications  to  MODFLOW  which  runs  the  forward 
simulation  for  the  groundwater  flow  problem.  That  is,  when  a  hydrogeologic  system  is  discretized  into 
a  model  structure  with  hydraulic  parameters  and  boundary  conditions,  MODFLOW  solves  for  head  or 
drawdown  and  flux  distributions  in  time  and  space.  If  the  model  structure  is  reasonable,  values  of 
parameters  are  accurate,  and  the  initial  and  boundary  conditions  are  defined  correctly,  one  may  expect 
that  the  calculated  solutions  (head  and  flux)  will  match  or  closely  approximate  the  observed  head  and 
flux.  However,  it  is  exceedingly  rare  to  have  enough  data  to  determine  “true”  values  of  aquifer 
parameters  for  a  hydrogeologic  system,  and  boundary  conditions  are  even  more  difficult  to  determine. 
As  a  result,  calculated  results  usually  deviate  from  the  observed  data.  Satisfactory  modeling  results 
may  often  be  obtained  only  after  many  trial-and-error  changes  to  system  parameters  and  boundary 
conditions.  The  changes  usually  rely  on  the  modeler’s  knowledge  of  hydrogeology,  degree  of 
familiarity  with  the  simulated  area,  and  professional  judgment.  The  decision  about  when  a  satisfactory, 
or  the  best,  model  has  been  achieved  largely  depends  on  the  modeler’s  subjectivity.  The  trial-and-error 
process  commonly  is  time-consuming  as  well. 

Alternatively,  parameter  estimation  may  be  achieved  by  using  nonlinear  regression.  Various 
methods  for  groundwater  parameter  estimation  have  been  developed  (e.g.,  Neuman  and  Yakowitz, 
1979;  Yeh,  1986;  Cooley  and  Naff,  1990;  Hill,  1992;  Sun,  1996).  We  have  been  working  with 
MODFLOWP  (Hill,  1 992)  for  estimation  of  parameters  in  a  three-dimensional  heterogeneous  aquifer 
system  by  modeling  aquifer  test  results  using  MODFLOW  for  the  forward  modeling.  In  addition  to 
the  utility  programs  for  MODFLOW,  we  have  developed  several  modifications  to  the  MODFLOWP 
code  to  enhance  usability.  In  particular,  we  have  added  changes  to  the  code  which: 

1 )  simplify  input  by  making  all  input  to  the  program  free  format; 

2)  allow  the  user  to  easily  change  the  number  and  combination  of  parameters  for  estimation  in 

a  given  run  by  modifying  one  line  rather  than  regenerating  the  whole  .PAR  file; 
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3)  automatically  calculate  ROFFs,  COFFs,  and  TOFFs  for  head  observations; 

4)  print  a  message  to  the  screen  when  estimated  parameter  values  are  outside  the  upper  and 

lower  limits;  and 

5)  allow  the  user  to  monitor  run  time. 

Each  of  the  MODFLOWP  utilities  is  discussed  below  with  an  example  of  its  usage. 

3.1.  Free  Format  for  Input  of  Data  to  the  Parameter  Estimation  Package 

With  previous  versions  of  MODFLOWP,  all  data  input  had  to  be  entered  according  to  highly 
structured  formats.  Data  input  to  the  .PAR  file  is  complicated,  and  the  formats  vary  between  data  sets 
in  the  file.  To  simplify  data  input  to  the  .PAR  file  and  reduce  errors  due  to  format  problems, 
MODFLOWP  has  been  modified  to  read  the  .PAR  file  as  free  format.  Data  sets  still  must  follow  the 
same  line  structure  (i.e.,  same  data  elements  on  a  given  line,  same  sequence  of  lines  and  data  sets)  and 
use  the  same  data  types  (i.e.,  integer,  real)  as  before.  However,  data  elements  on  a  given  line  can  be 
given  in  free  format,  and  so  they  need  only  be  separated  by  one  or  more  spaces. 

Structured  format  entry  interprets  blank  spaces  as  zeros  leading  to  a  practice  of  not  providing 
data  for  input  values  that  are  intended  to  have  a  value  of  zero.  Free  format  requires  zero  values  to  be 
specified.  Hence,  all  data  values,  including  zero,  must  be  supplied  when  using  free  format  input. 

3.2.  Enable  Any  Number  and  Combination  of  Parameters  for  Estimation 

MODFLOWP  input  file  format  requires  that  parameters  to  be  estimated  should  precede  those 
that  will  remain  unchanged.  In  Data  Set  9  all  positive  group  numbers  need  to  precede  any  negative 
numbers  (Hill,  1992).  If  the  MODFLOWP  user  wishes  to  evaluate  different  combinations  of 
parameters,  the  .PAR  file  has  to  be  reorganized.  Rebuilding  the  .PAR  file  can  be  time-consuming  and 
may  lead  to  input  errors.  The  modified  version  of  MODFLOWP  presented  here  allows  the  user  to 
choose  any  new  or  former  combination  of  initially  identified  parameters  in  Data  Set  9  for  estimation 
without  reorganizing  the  .PAR  file.  The  user  simply  identifies  in  Data  Set  9  those  parameters  that  will 
not  be  evaluated  by  assigning  negative  group  numbers  to  them.  The  negative  numbers  can  be  in  any 
sequence. 

3.3.  Calculate  and  Enter  Offsets  inside  MODFLOWP  _  r  -  -  -- 

MODFLOWP,  as  a  numerical  model,  operates  by  discretizing  time  into  time  steps  and  space 
into  grid  cells  with  nodes  in  the  centers  of  cells.  Calculations  of  model  results  are  made  at  time  steps 
and  nodes  which  generally  do  not  correspond  exactly  with  the  time  of  observations  or  the  locations  of 
well  screens.  Such  differences  may  not  be  significant  in  many  modeling  tasks,  but  they  could  add 
model  error  to  the  parameter  estimation  process  and  they  may  be  significant  for  some  modeling 
scenarios,  especially  for  modeling  of  pumping  test  results  at  closely  spaced  partially  penetrating  wells 
in  heterogeneous  media. 

MODFLOWP  allows  the  user  to  specify  offsets  between  locations  of  observation  points  and 
nodes  setting  the  ROFF  and  COFF  data  values  in  the  .PAR  file,  and  between  observation  times  and 
model  time  steps  by  setting  TOFF  values.  The  user  must  calculate  each  ROFF,  COFF  and  TOFF 
outside  of  the  code.  These  offsets  must  be  recalculated  and  reentered  each  time  the  time  or  grid 
discretization  are  altered  in  subsequent  model  runs.  Both  activities  can  be  time-consuming  and  error- 

*  If  an  F77  compiler  is  used, 
character  input  may  need  to  be 
bracketed  by  apostrophes.  For  23 

example  in  Example  File  16, 
line  10,  Q  must  be  input  as  ’ Q ’ . 


prone.  We  have  incorporated  utilities  into  MODFLOWP  which  automatically  calculate  these  offsets 
eliminating  the  need  to  enter  them  into  the  .PAR  file. 
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Figure  4.  ROFF  and  COFF  defined  by  absolute  coordinates  (x0,  y0).  Origin  of  x-y 
portion  of  coordinate  system  is  in  the  upper  left-hand  corner  of  the  model  domain 
(the  GMS  default  coordinate  system). 

3.3.1.  Calculate  and  Enter  ROFFs  and  COFFs 

ROFFs  and  COFFs  represent  the  offsets  of  an  observation  well  from  the  center  of  the  cell  in 
row  and  column  directions,  respectively.  They  are  used  to  calculate  the  head  (or  drawdown)  for  the 
observation  well  by  interpolation  of  the  calculated  heads  (or  drawdown)  of  nearby  cells.  In 
MODFLOWP,  ROFFs  and  COFFs  have  to  be  provided  in  Data  Set  6  of  the  .PAR  file.  If  grid 
discretization  or  locations  of  observation  points  are  changed,  ROFFs  and  COFFs  have  to  be  changed 
also. 

Although  the  calculation  of  ROFFs  and  COFFs  is  not  difficult  (see  MODFLOWP  user’s  guide 
[Hill,  1992]),  it  is  simpler  to  directly  input  the  coordinates  (Xq,  y0)  of  an  observation  well  (Figure  4) 
and  have  the  program  calculate  ROFFs  and  COFFs  internally.  That  is,  (x0,  y0),  the  plane  coordinates 
of  an  observation  well  are  input,  rather  than  ROFF  and  COFF  and  the  row  and  column  numbers  of  the 
well.  There  are  at  least  two  advantages  with  this  input  approach:  ( 1 )  x,  y  coordinate  locations  are 
accurate  and  intuitive,  and  (2)  ROFFs  and  COFFs  will  be  automatically  recalculated  inside  the 
program  whenever  the  spatial  discretization  is  changed. 

3.3.2.  Calculate  and  Enter  TOFFs 

TOFFs  (time  offsets)  in  MODFLOWP  are  used  to  calculate  differences  between  actual 
observation  times  and  the  time  steps  used  in  the  simulation.  Standard  MODFLOWP  requires  the  user 
to  input  TOFFs  and  the  corresponding  time  step  in  Data  Set  6.  Calculation  and  input  for  TOFFs  are 
not  difficult  but  can  be  time-consuming.  Whenever  temporal  discretization  is  changed,  the  time  steps 
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and  TOFFs  for  Data  Set  6  and  Data  Set  6C  have  to  be  changed. 

Our  modifications  have  made  this  task  very  simple:  just  input  observation  time  and  observed 
head  H;,  then  the  corresponding  time  steps  and  TOFFs  are  internally  calculated  within  MODFLOWP. 
Advantages  for  direct  input  of  t;  and  H,  are:  (1)  users  are  most  familiar  with  the  measured  head  at  a 
specific  time  making  input  errors  easier  to  spot;  (2)  records  of  field  measurements  are  usually  time  and 
head  (or  drawdown);  thus  the  field  records  can  be  used  directly  with  our  modified  input  format  and 
compared  directly  with  the  input  file. 

3.3.3.  Input  Instructions  for  Data  Set  6 

Corresponding  with  the  modifications  described  above,  the  inputs  for  Data  Set  6  and  Data  Set 
6C  need  to  be  changed  as  follows: 

Data  Set  6:  DID,  NDER,  X  ,Y,  OBST,  HOBS,  WT1,  WT2,  IWT,  IPUT 

•  DID  and  NDER  are  the  same  as  defined  in  the  MODFLOWP  user’s  guide. 

•  X  and  Y  are  the  coordinates  of  an  observation  well. 

•  OBST  is  the  time  an  observation  was  made.  OBST  >  0  is  the  observation  time  if  there  is  only 
one  observation  at  this  well.  Otherwise  OBST  <  0  and  the  absolute  value  of  OBST  equals  the 
number  of  observations  at  this  well,  and  OBST  must  be  entered  as  a  real  number.  For 
example,  if  three  observations  are  made,  OBST  =  -3.0. 

•  Definitions  for  HOBS,  WT 1 ,  WT2  and  IWT  are  the  same  as  in  the  MODFLOWP  user’s  guide. 

•  IPUT  is  a  new  flag  defined  in  this  report  indicating  whether  the  same  weights  (IPUT  =  0)  will 
be  used  for  all  observations  at  this  well,  or  if  different  weights  (IPUT  >  0)  will  be  specified. 


Data  Set  6C:  DID,  OTIME,  HOBS,  WT1,  WT2,  IWT 

•  DID  is  the  same  as  defined  in  the  MODFLOWP  user’s  guide. 

•  OTIME  is  the  observation  time.  This  replaces  time  step  and  time  offset  in  Data  Set  6C 

•  Definitions  for  HOBS,  WT  1 ,  WT2  and  IWT  are  the  same  as  in  the  MODFLOWP  guide. 

3.4.  Announce  Parameter  Values  out  of  Upper  and  Lower  Limit  Bounds 

In  MODIUOWP,  parameter  bounds  are  only  used  for  reference  in  the  output  file  and  are  not 
used  in  the  regression  solution.  Estimated  parameters  may  be  outside  the  bounds  of  expected, 
reasonable,  bounds.  Lor  example,  during  a  MODFLOWP  ran  specific  yield,  Sy,  may  be  estimated  to 
be  larger  than  0.5;  specific  storage,  Ss,  may  become  zero;  K  for  a  sand  layer  may  have  an  estimated 
value  that  is  lower  than  that  of  a  shale. 

To  identify  parameter  estimation  runs  that  are  in  progress,  but  not  likely  to  result  in  a 
meaningful  set  of  parameter  estimates,  a  message  will  be  sent  to  standard  output  (screen)  if  any 
evaluated  parameter  is  set  beyond  one  of  the  bounds  specified  in  Data  Sets  8 A  and  8B.  The  message 
identifies  the  parameter,  the  parameter  value,  and  the  parameter  limit  violated.  This  allows  the 
modeler  to  terminate  the  MODFLOWP  run,  which  may  take  hours  or  days  to  complete,  when  it  is 
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obvious  the  program  is  not  converging  on  a  reasonable  solution.  A  message  will  also  be  printed  if  the 
value  of  a  parameter  changes  sign  during  an  estimation  run. 

3.5.  Monitor  Run  Time 

As  mentioned  above,  a  given  parameter  estimation  run  may  take  a  long  time  because  of  the 
number  of  parameters  being  estimated,  the  complexity  of  the  hydrological  processes  being  modeled, 
and  the  nonlinear  nature  of  the  problem.  For  the  modeler,  feedback  on  the  amount  of  CPU  time  used 
for  a  given  iteration  is  sometimes  a  useful  indication  of  whether  the  investigated  problem  is 
“moderately”  or  “highly”  nonlinear,  because  a  highly  nonlinear  problem  usually  takes  a  longer  time  to 
converge.  The  modified  program  prints  to  file  and  writes  out  to  the  screen  the  elapsed  run  time  since 
the  starting  date  and  time  for  each  iteration. 


Monitoring  the  run  time  requires  the  use  of  non-standard  FORTRAN  routines.  We  have 
supplied  seven  versions  of  the  run  time  routines.  The  following  table  lists  file  names  to  be  used  for 
specific  compilers: 


time_f90.for 

time_g77.for 

time_ms.for 

time_duf.for 

time_sun.for 

time_sgi.for 

time_dum.for 


FORTRAN  90 

GNU  g77  FORTRAN 

Microsoft  Power  Station  FORTRAN 

Digital  FORTRAN  for  UNIX 

Sun  FORTRAN 

SGI  FORTRAN 

Disables  time  calls  for  use  with  any  compiler 


The  list  of  compiler  coverage  is  far  from  exhaustive.  However,  many  compilers  adhere  to  non- 
FORTRAN  standards.  Hence,  the  coding  is  more  generally  applicable  than  listed  here.  The  routine 
time_f90.for  uses  the  FORTRAN  90  standard  time  calls  which  are  supported  by  a  large  number  of 
FORTRAN  77  and  FORTRAN  90  compilers.  Using  time_dum.for  removes  time  monitoring  from  the 
program  without  requiring  modification  of  the  source  code.  The  file  time_dum.for  can  be  used  with 
any  compiler  because  it  simply  returns  to  the  call  statement  without  executing  the  timing  routines. 

3.6.  Example  .PAR  and  Output  Files 

The  following  example  is  taken  from  Test  Example  1  of  the  MODFLOWP  user’s  guide  (Hill, 
1992),  but  is  changed  for  use  with  the  modified  version  of  MODFLOWP  presented  with  this  report. 
Consider  the  sequence  of  modeling  runs  where  initially  9  parameters  are  estimated  using  the  observed 
drawdown  from  several  wells.  In  our  test  of  the  modified  MODFLOWP  program,  we  start  with  the 
same  input  sequence  of  9  parameters  but  now  we  just  want  to  estimate  7  of  them:  in  sequence 
positions  2,  4,  5,  6,  7,  8,  and  9.  To  evaluate  parameters  in  sequence  positions  2,  4,  5,  6,  7,  8,  and  9 
only,  the  user  changes  the  sign  in  Data  Set  9  for  the  group  numbers  of  those  parameters  that  will  NOT 
be  estimated:  -3,  4,  -5,  6,  3,  4,  7,  7,  1  (Example  File  16).  In  Example  File  16,  Comments  lines  starting 
with  @@@  have  been  added  to  the  .PAR  file  listing  to  highlight  the  changes  made  to  the  .PAR  file  for 
Test  Example  1.  In  addition  to  the  modification  of  Data  Set  9,  observation  well  drawdown  data  and 
the  pumping  well  drawdown  data.  Data  Set  6C,  has  been  input  with  free  format. 
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The  affected  portions  of  the  corresponding  output  file  are  given  in  Example  File  17. 


Example  File  16.  Annotated  (@@@)  portions  of  the  MODFLOWP  .PAR  file. 


TWO-LAYER  EXAMPLE  -  TRANSIENT  LINE  1 

Modified  Test  Case  1  of  Hill  (1992 .Appendix  A)  LINE  2 


9 

5 

56 

4 

2  3 

0  1 

2 

1 

1 

32 

0 

0 

1 

18  3 

0  0 

1 

-1 

1 

0 

0 

40 

33 

0 

0 

0  0 

0  0 

@@@ 

All 

zero 

values 

must  be 

included 

1 

0 

0 

0 

30 

Q  0 

-2 

l 

15 

-1  0 

1. 

0 

1 

9 

10 

1.0 

2 

9 

10 

1.0 

SI  1 

-1 

0 

0 

0  0 

@@@ 

All 

zero 

values 

must  be 

included 

LINE  3 
LINE  4 
LINE  5 
LINE  6 

0  0  0  0  0  LINE  7 


LINE  8 
DATA  SET  1A 
DATA  SET  2 

DATA  SET  2A 

DATA  SET  2 


1-0  1  500.  2500.  -3.0  0.00  1.0025  0.0025  0  0 

@  @  @  Input  x  and  y  coordinates.  Number  of  observations  is  real  value 


DATA  SET  6 


.  „  *  .  DATA  SET  6B 

1.0  0.00  101.804  DATA  SET  6C 

@  @  @  Input  observation  time  not  time  step  and  offset 
@  @  @  Observation  weighting  taken  from  Data  Set  6 

1.1  87162.0  101.775 

1.12  24439068.0  101.675 


6-0  2  3500.0  3500.0  -3.0  0.0  0.000  0.000 

@@@  transient  observation  weight  option  set  to  1 

2 

6.0  0.00  126.537  1.0025  0.0025  0 

@@@  Full  control  of  weighting  for  each  observation 


1  DATA  SET  6 


DATA  SET  6B 
DATA  SET  6C 


6.1 

87162.0 

126.542  1. 

0025  0.0025 

0 

6.12 

24439068.0 

112.172  1 

.0025  0.0025 

0 

-1.074 

1 . 38E-3 

. 426E-03 

1.2E-3 

t.E-7 

4.0E-5 

2.E-8 

1.0E-8 

-0.8 

1.3E-2 

3.00E-3 

1.2E-2 

l.E-6 

4 . OE-4 

2.E-7 

1.0E-7 

-1.4 

1.3E-4 

3.00E-5 

1.2E-4 

l.E-8 

4.0E-6 

2.E-9 

1.0E-9 

-3 

4 

-5 

6 

3 

7 

7 

1 

@  @  @  inactive  parameters  may  be  interspersed  with  active  parameters 


2  .  E-4 

DATA  SET  8 
2  .  E-3 

DATA  SET  8A 
2  .  E-5 

DATA  SET  8B 

4 

DATA  SET  9 


Example  File  17.  Annotated  (  @@@  )  portions  of  the  MODFLOWP  output  file. 

1  MODFLOW 

U.S.  GEOLOGICAL  SURVEY  MODULAR  FINITE-DIFFERENCE  GROUND-WATER  FLOW  MODEL  -  MODFLOWP 

Bsu  August,  1998  Modified  version  @@@  Announcement  of  special  version 

See  Boise  State  University  Technical  Report  BSU  CGISS  97-02 

OBSERVED  HEAD  DATA 
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@@@  Observation  time  or  number  of  transient  observations  as  item  6. 


OBS#  ID  LAYER  ROW  COLUMN  OBS .  TIME 
1  1.0  13  1  -3.000 

TRANSIENT  DATA  AT  THIS  LOCATION,  ITT  = 

1  i.o  .ooo  o 

2  1.1  87162.000  1 

3  1.12  24439070.000  12 


TIME  STEP  ROW/COLUMN/TIME  OFFSETS  OBS.  HEAD  STATISTIC 
-3  .000  .000  .000  .000  1.00250 

2 


.000 

.000 

.000 


101.80 

-.28999E-01 

-.12900 


1.0025 
. 25000E-02 
. 25000E-02 


@  @  @  Observation  time  and  time  step  and  offset  listed . 


PARAMETER  INFORMATION: 


ID 

t  EST 

INITIAL  VALUE 

LN 

UPPER  VALUE 

LOWER  VAI 

SI 

. 13  80E-02 

1 

. 1300E-01 

. 1300E-03 

KST 

. 1200E-02 

1 

. 1200E-01 

. 1200E-03 

KV 

. 1000E-06 

1 

. 1000E-05 

. 1000E-07 

SI 

. 2000E-03 

1 

. 2000E-02 

. 2000E-04 

T 

. 4000E-04 

1 

. 4000E-03 

. 4000E-05 

RCH 

. 2000E-07 

0 

. 2000E-06 

. 2000E-08 

RCH 

-1000E-07 

0 

. 1000E-06 

. 1000E-08 

Q 

-1.074 

0 

-.8000 

■1.400 

T 

. 4260E-03 

0 

. 3000E-02 

.  3000E-04 

CONVERGENCE  CRITERIA  GROUP#  WEIGHT 


@  @  @  inactive  parameters  listed  last . 


.  725E-02 
. 833E-02 
100. 

. 500E-01 

.250 

500. 

. 100E+04 

.000 

.000 


.000 
.  000 
.000 
.000 
.000 
.000 
.000 
.000 
.000 


4 

6 

3 

4 
7 
7 
1 

-3 

-5 


SCALED  SENSITIVITIES  (SCALED  BY  B*(WT**.5)) 


PARAMETER  #  : 
PARAMETER  ID 
OBS#  &  ID 
1  1.0 
2  1.1 


2 

SI 


4 

KST 


5 

KV 


.000  1.51  - . 445E-07 

. 123E-01  - . 3 97E-02  -.113E-02 


COMPOSITE  SCALED  SENSITIVITIES: 

((SUM  OF  THE  SQUARED  VALUES) /ND) **. 5 
235.  25.7  123. 


6 

SI 


.000 

• . 3  59E-02 


57.7 


7 

T 

-.228E-04 
. 621E-02 


498. 


8 

RCH 

.150 

.  800E-06 


11.9 


9 

RCH 

. 749E-01 
. 400E-06 


11.0 


Run  time  [min] : 


.1500  since  8/12,  11:24 


ITERATION  NO.  =  1 

VALUES  FROM  LEAST-SQUARES  REGRESSION  PROCEDURE  : 


DET  OF  SCALED  LEAST-SQUARES  MATRIX  =  .25575E-02 

MARQUARDT  PARAMETER  (AMP) - =  .00000 

FACTOR  FOR  SCALING  PAR.  CHANGE  (AP)=  1.0000 

MAX.  FRACTIONAL  PAR.  CHANGE  (DMX)  =  .76297 

MAX.  FRAC.  CHANGE  OCCURRED  FOR  PAR.#  5 

UPDATED  ESTIMATES  OF  REGRESSION  PARAMETERS  : 

SI  KST  KV  SI 

PAR.#  2  4  5  6 

0  . 11753E-02  . 1133  IE-02  .17630E-06  .10834E-03 


@@@  Run  time  of  first  iteration 


T 

RCH 

RCH 

7 

8 

9 

49218E-04 

.  11507E-07 

.  15596E-07 

Run  time  [min]:  .2167  since  8/12. 11:24 

ITERATION  NO.  =  2 


PARAMETER  SUMMARY 


PARAMETER  #  : 

2 

4 

5 

PARAMETER  ID  : 

SI 

KST 

KV 

FINAL  VALUES 

- .  677E+01 

- . 682E+01 

- . 153E+02 

EXPONENTIAL 

OF  LN  PARAMETERS 

(0.0  FOR  UNTRANSFORMED  PARAMETERS) 

. 115E-02 

. 109E-02 

. 222E-06 

STD .  DEV . 

. 387E-01 

. 7  97E-01 

. 7  82E-01 

COEF.  OF  VAR. 

. 572E-02 

. 117E-01 

. 510E-02 

UPPER  LIMIT,  LINEAR  95-PERCENT  CONF .  INT. 

- . 669E+01  - . 666E+01  -.152E+02  - 
EXPONENTIAL  OF  LN  PARAMETERS 


@@@  Cumulative  run  time  at  each  iteration 


6 

7 

8 

9 

SI 

T 

RCH 

RCH 

988E+01 

- . 994E+01 

. 109E-07 

.  160E-07 

512E-04 

. 481E-04 

. 000E+00 

. OOOE+OO 

706E+00 

. 513E-02 

. 128E-08 

. 680E-09 

715E-01 

. 516E-03 

. 118E+00 

. 426E-01 

847E+01 

- . 993E+01 

. 134E-07 

. 173E-07 
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(0.0  FOR  UNTRANSFORMED  PARAMETERS) 

. 125E-02  . 128E-02  .259E-06  .210E-03  .486E-04  .OOOE+OO  .000E+00 

LOWER  LIMIT.  LINEAR  95-PERCENT  CONF.  INT. 

- . 684E+01  - . 698E+01  -.155E+02  -.113E+02  -.995E+01  .830E-08  .146E-07 

EXPONENTIAL  OF  LN  PARAMETERS 
(0.0  FOR  UNTRANSFORMED  PARAMETERS) 

. 107E-02  . 932E-03  .190E-06  .125E-04  .476E-04  .OOOE+OO  .000E+00 

Runtime  [min]:  .4500  since  8/12,  11:25  @  @  @  Total  mn  time  USed. 

4.  OBTAINING  SOURCE  CODE 

Source  code  of  our  modifications  is  available  under  the  listing  for  this  report  at 
http://kihei.idbsu.edu/cgiss_pub.html.  Compacted  archives  in  both  WINZIP®  and  tar  formats  are 
available  which  include  a  modified  source  code  for  both  MODFLOW  and  MODFLOWP,  example 
input  and  output  files,  and  UNIX  ‘diff  files  which  highlight  the  changes  made  to  the  MODFLOW  and 
MODFLOWP  programs.  To  obtain  a  copy  of  these  archives  connect  to 

http://kihei.idbsu.edu/cgiss_pub.html  with  a  web-browser,  follow  the  Technical  Report  Publications 
link,  and  then  scroll  to  Fluang  et  al..  The  report  listing  has  been  set  up  as  a  link  to  downloading 
instructions  for  the  source  code  archive  files. 

5.  SUMMARY 

This  report  serves  as  documentation  and  user's  guide  for  ( 1 )  a  number  of  stand-alone  utilities 
to  assist  with  input  to  MODFLOW  and  (2)  code  modifications  internal  to  MODFLOW  and 
MODFLOWP.  These  were  developed  to  improve  efficiency  of  use  of  MODFLOW  and  MODFLOWP 
for  modeling  groundwater  flow  in  three-dimensional  heterogeneous  systems.  Several  utilities  to  assist 
with  input  to  MODFLOW  are  specifically  designed  for  use  with  the  pre-  and  post-processor  package: 
Groundwater  Modeling  System  (GMS)  and  are  identified  as  such  in  this  report.  New  features  for 
MODFLOW  include  routines  for: 

1 )  assignment  of  material  blocks  and  material  types  to  these  blocks  within  a  given  three- 

dimensional  grid  mesh; 

2)  distribution  of  pumping  rate  for  individual  model  layers  for  a  well  that  draws  from  or  injects 

into  multiple  layers; 

3)  calculation  of  weighted  drawdown  for  an  observation  well-  screened-over  multiple  layers; 

and 

4)  customizing  the  drawdown  output  file  to  contain  data  for  a  specified  number  of  observation 

wells  (rather  than  for  all  nodes  individually). 

For  MODFLOWP,  the  program  has  been  modified  to: 

1)  simplify  input  by  making  all  input  to  the  program  free  format; 

2)  allow  the  user  to  easily  change  the  number  and  combination  of  parameters  for  estimation  in 

a  given  run  by  modifying  one  line  rather  than  regenerating  the  whole  .PAR  file; 

3)  automatically  calculate  ROFFs,  COFFs,  and  TOFFs; 

4) announce  parameter  estimates  outside  the  upper  or  lower  limits  during  program  execution; 

and 

5)  allow  the  user  to  monitor  run  time. 
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